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We consider deterministic and randomized quantum algorithms simulating e~~' Ht by a product of 
unitary operators e " iAjtj , j = 1, . . . ,N, where Aj € {Hi, H m }, H = YT=i H * and *i > for 
every j. Randomized algorithms are algorithms approximating the final state of the system by a 
mixed quantum state. First, we provide a scheme to bound the trace distance of the final quantum 
states of randomized algorithms. Then, we show some randomized algorithms, which have the 
same efficiency as certain deterministic algorithms, but are less complicated than their opponentes. 
Moreover, we prove that both deterministic and randomized algorithms simulating e~ lHt with error 
£ at least have f2(t 3/,2 £ -1 / 2 ) exponentials. 
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I. INTRODUCTION 



While the computational cost of simulating many particle quantum systems using classical computers grows expo- 
nentially with the number of particles, quantum computers nonetheless have the potential to carry out the simulation 
efficiently [6( . This property, pointed out by Feynman, is one of the founding ideas of the field of qua ntum computation. 
The simulation problem is also related to quantum walks and adiabatic optimization 0, [l(| [ll|, [l2|, HH • 

A variety of quantum algorithms have been proposed to predict and simulate the behavior of different physical and 
chemical systems. Of particular interest are splitting methods that simulate the unitary evolution e~ lHt , where H is the 
system Hamiltonian, by a product of unitary operators of the form e~ tAjtj , j = 1, . . . , N, where Aj £ {Hi, . . . , H m }, 
i ^ , '' H = X)"=i H an d assuming the Hi do not commute. 

A recent paper [l| shows that high order splitting methods [3, 0] can be used to derive bounds for N that are 
7-H ■ asymptotically tight. This work also provides algorithms that achieve the upper bounds. However, the derived 
algorithms require some of the tj to be negative, which may limit their application. For instance, the algorithms 
cannot be used for the simulation of diffusion operators, because there exists no inverse exponential diffusion operator, 
as noted by Suzuki j3] who proposed the high order splitting methods. The reason is that for spiting methods with 
order of convergence greater than or equal to three, some of the {tj} must be negative 

In this paper, we consider deterministic and randomized quantum algorithms simulating e~ lHt using only positive 
{tj}. By randomized algorithms we mean algorithms approximating the final state of the system by a mixed quantum 
state. We show that: 



1. The increase in the trace distance of the quantum states in a randomized algorithm is bounded from above by 

2 || E(U^)-U Q || +E(\\ U u -U f), 

where Uo is the unitary evolution being simulated, U u denotes the randomized algorithm and E(-) is the 
expectation. 

2. Deterministic and randomized algorithms simulating e~ lHt by approximating it by Ilj=i e~ lAjtj with error e 
must satisfy 

3. The optimal deterministic algorithm is based on the Baker-Campbell-Housdorf formula Q. 

4. An optimal randomized algorithm is obtained by a direct application of the Trotter formula 0- 



'Electronic address: cz2165@columbia.edu 



2 



II. RANDOMIZED ALGORITHMS FOR QUANTUM SIMULATION 

Let us now state the problem in more details, then discuss the algorithms and their performance. A quantum 
system evolves according to the Schrddinger equation 

ij t m))=H\m), (i) 

where H is the system Hamiltonian. For a time-independent H, the solution of the Schrddinger is 

|^))=e- lHt |Vo), (2) 
where \tpo) is the initial state at t = 0. Here we assume that H is the sum of local Hamiltonians, i.e., 



H = Y,H k , (3) 



fc=i 



and all the Hk are such that e tHhT can be implemented efficiently for any r > 0. Therefore, we will be using a 
product of the unitary operators U = Y\.j=i e~ lAjtj , where Aj 6 {Hi, . . . , H m }, tj > 0, to simulate Uo — e~ lHt . 
However, since the Hk do not commute in general, it introduces an error in the simulation. We measure this error 
using the trace distance, as in Out goal is to obtain tight bounds for N for algorithms achieving accuracy e in 
the simulation and to show optimal algorithms. 

Variations of this problem that do not set restrictions on tj have been extensively studied in the literature, see, 
e.g., 0,11,0. As far as we know only deterministic algorithms have been considered. In this paper, we propose a 
randomized model for simulating quantum systems, which simplifies the design of algorithms without compromising 
their efficiency. 

In the randomized model, the sequence of unitary operators is selected randomly according to a certain probability 
distribution. The distribution can be realized either by "coin-flips" or by "control qubits" , which requires some 
ancillary qubits. As a result, the algorithm is a product of a random sequence of unitary operators U w = Il^i e~ iA ^> u ' t i 
selected with probability p^ . Hence, the final state of the quantum algorithm is the mixed state 

p = J2pM^o)(MuI (4) 

For more general cases, where the input state of the simulation is not exactly |^o)> but a different mixed state po, the 
final state is p = J2^P^ u ^PoUl- 

In order to analyze the efficiency of randomized algorithms, we show an upper bound for the trace distance between 
the desired final state and the one computed by a randomized algorithm. 

Lemma 1. Let Uo be the unitary evolution being simulated by a set of random unitary evolutions as we described 
above. Then the trace distance between a = J7o|i/'o)('0o|^o an( ^ P * s bounded from above by 

D(po,\iM(M)+ 2 \\ 52puUu,-u ii +^>c II u u -u 1| 2 

uj u (5) 

=D(p , IVoXVoI) + 2 || E(UJ) - Uo || +E(\\ U u - U f), 

where D(-) denotes the trace distance and E(-) denotes the expectation. 

Proof. First, we calculate the difference of the output states pi and p\, which is 

Y,P^PoUl-U Q \^o)(MU^ 
= J2(Uo + U„- Uo)po(U + U„- Uo) ] - WoXV-ol^ 

UJ 

= Y,P^ - U °)P° U l +J2p^Po(U^ - Uo) +Y,P"( U " - Uo)po(U u - [/ ) t +5>u,t/ />ot4 - WoXV-olE/o* 

UJ UJ UJ UJ 

=(52p<*U u - Uo) P oUl + UoPoC^Tp^ - - U )po(U ul - Uo) + U ( P o ~ |Vo><Vo|)f4- 

UJ UJ UJ 

(6) 
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Hence, 

D(pi,pi) =Tr|pi- Pl | 

<Tr|(^ Pw [/ w - Uo) P oUl\ + Tr\Uopo(52p<*U u - t/o) f | + ^p w Tr|([/ w - E/b)po(£k - U )\ + Tr|p - \ih)(M 

<2 || ^p u U u -Ua\\+J2P" II U " - ^ II 2 + Tr l/°o - IVoXlMI 
=D(po, |^o)(^o|) + 2 || £(*7 W ) - C/ II Ua, - U || 2 ). 

□ 

Let AD — D(pi,p[) — D(p 0} V'oKV'ol), which is the augment of trace distance. From the above lemma, we know 
in the simulation, AD is bounded from above by || — Uo | and E(\\ — Uo || 2 ). Moreover, it is easy to check 

that AD = 6(|| E(U U ) - U ||) for certain p , as well as AD = Q(E(\\ U u - U || 2 )) for some other p . Therefore, the 
lower bound is also tight asymptotically, i.e., 

AD = 6(2 || E(U U ) - U || +E(\\ U w - U a || 2 )). 

For the convenience of the reader, below we give two examples of randomized algorithms and we use the lemma 
above to analyze their cost. It turns out that the second algorithm is optimal. 

• Algorithm 1 

Divide the total evolution time t into equal K small segments of size At. 
Let po = \ipo)(ipo\ be the input to the first stage of the algorithm. 

Consider the fc-th stage of the algorithm where the input is Pk-i, k — 1, . . . ,mK. The algorithm chooses 
uniformly and independently at random operators from {e~ zHlAt , . . . , e -*#™ A *}. Hence, the output of stage k 
is 

m 



The final result of the algorithm is p m x and is used to approximate a. 

Due to Lemma 1, the error of this algorithm for simulating e ~ lHAt by m consecutive stages (i.e., by the stages 
km + 1 and (k + l)m, for any k = Q,...,K—l) is bounded by two elements, || E(U U ) — Uo || and E(\\ U u — Uo || 2 ), 
where U u is the product of the sequence m operators. Since the selection in each stage is independent and uniform, 

1 m m 

E{U„) = (- H~ iH ™ M ) m =I-iJ2 HjAt + 0(At 2 ). (8) 

m 3 = 1 3 = 1 

Hence, || E(U„) - U Q \\= 0{At 2 ). Furthermore, for any u, U u = I + O(At), then E(\\ U u - U || 2 ) = 0{At 2 ). 
Therefore, the error in each m consecutive stages is 0{At 2 ). Thus, the total error is e = 0{KAt 2 ) and the total 
number of exponentials used is TV = mK = 0(t 2 je). 

We remark that this is equal modulo a constant to the cost of the deterministic algorithm that is based on a direct 
application of the Trotter formula, i.e., the one that uses 

m 

If e- iH ' At 

3 = 1 

to simulate e~ lHAt . However, Algorithm 1 has certain advantages over this deterministic algorithm. In the deter- 
ministic algorithm, in order to simulate e ~ lHAt , we need to store the current index j of e ~ lH ' At , for j = 1, • • • , m. 
However, in Algorithm 1, each stage is independent and the algorithm is "memoryless" . 

• Algorithm 2 

Divide the total evolution time t into equal K small segments of size At. 
Let po — IV'oXV'ol be the input to the first stage of the algorithm. 
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Consider the fc-th stage of the algorithm where the input is pk-i, k — 1, . . . , K. The algorithm select an operator 
uniformly and independently at random from the set of operators 

m 

{J^J e ~ lH "U) At : (j varies over all permutations on m symbols }. 
Then, the output of the fc-th stage is 



p* = E i(Q - ' v ^ i ll <""" u • 1 



The final result of the algorithm is and is used to approximate a. 



Let 



e -iH a (j)At 



= II^ - ^-Ci) A< - o^0) A * 2 + °( A * 3 )) 

= 7 - * E ff *Ci) A * - 2 E ^0) A * 2 -EE H„ {j) H a (k)At 2 + 0(At 3 ) 

j=l j=l j<kj<k 

TTL 7YI 

= I- l J2H 3 At-lj2 H l At2 ~ J E H„ U) H a (k)At 2 + 0(At 3 ). 



2 ^ 3 2 

j=l j=l j<k 

In each stage, the simulating operator XJ U could be U a with probability 1/m!, for each a in the permutations on 
m symbols. Since, for each a, \\ U a - U a \\= 0(At 2 ), so E(\\ U w - U || 2 ) = 0(At 4 ). Moreover, E(U U ) = I - 
iJ^LiHjAt - l^liHjfAt 2 + 0(At 3 ), hence || E{U W ) - U \\= 0(At 3 ). Due to Lemma 1, the error of this 
algorithm for simulating e - lHAt i n each stage is 0(At 3 ). Thus the total error e = 0(KAt 3 ). Hence, for a given e 
the value of K in Algorithm 2 is smaller than that in Algorithm 1. The total number of exponentials used is 
N = mK = 0(i 3 / 2 /£ 1/2 )- 

We remark that this is equal modulo a constant to the cost of a deterministic algorithm solving the problem. The 
difference is that the deterministic algorithm is more slightly more complicated than the one discussed in the previous 
item. It is based on the Baker-Campbell-Housdorf formula (Strang splitting) and uses 

m 1 



j=l j=m 



to simulate e lHAt . 



III. LOWER BOUNDS FOR RANDOMIZED ALGORITHMS 



In fact, Algorithm 2 is asymptotically optimal among all randomized algorithms simulating the evolution of the 
quantum system. Before proving the optimality of Algorithm 2, we start with a lemma. 

Lemma 2. For any < Xi < 1, i = 1, • • • , N, and x i = 2, let S be the sum of all elements in {xiXjXk : i < 

j < k,2\k-i,2\j-i}. Then S< 1/3. 

Proof. For N = 3, it is easy to check S < (2/3) 3 < 1/3. 
Assume that, for N < M, the conclusion holds. 

For the case N — M, the global minimum of S will be achieved at a local minimum or the border. Moreover, if the 
global minimum is obtained at the border, which means some Xi = or 1, it reduces to the case N < M. Then, we 
only consider its local minimum and assume Xi ^ for each i. 
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Let f(x) = S — A(^^ Xi — 2). In any local minimum, J^- = 0, for each Xi. Then 



dx. 



x 3 x k + Xj Xk + XjX k - A = 



i<j<k j<i<k j<k<i 

2\j-i 2\j-i 2\j-i 

2\k-i 2\k-i 2\k-i 



, and 

df 



dx, 



+ 2 



XjXk + ^2 x 3 xk + X] xjXk _ ^ = o- 

i+2<j<k j<i+2<k j<k<i+2 

2\j-i 2\j-i 2\j-i 

2\k-i 2\k-i 2\k-i 



Combined these two equations, we have 

Xi+i( x k + ^ Xk - ^ Xk - Xk) = 0. 

k>i+2 k<i-l k>i+3 k<i 

2\k-i 2\k-i 2\k-i 2 1 fc — i 

From the assumption, Xi+i ^ 0, we have 

h Xi-3 + Xi-i + X l+2 + X l+4 H = h Xi-2 + Xi + x i+3 + x l+5 H . 

Then, we consider — and d ^ +3 — 0, which can derive 

h Xi-l + X i+1 + X i+4 + X l+5 + • • • = h .Xj_2 + Xi + x i+3 + x i+5 -\ . 

Combine them together, we have x i+1 = x i+2 . Therefore, x 2 = x 3 = ■ ■ ■ = x N _i. Then, by considering = 
and Sr- = *S f 7 we have x\ = when TV is even; and x\ = x 2 = ■ ■ ■ = xn when N is odd. Since the first case 
is contradict to our assumption, we need only consider the case N is odd. Let TV = 2K + 1, then each term in S is 
(j^) 3 , and there are ^K(K + 1)(2K + 1) terms. Therefore, at the local minimum 

s = \ K ^ K + W + ^ JtT )3 = ^ - ^ < T 

Hence, for any N the conclusion holds. □ 

From the above Lemmas, we have the following two theorems. 

Theorem 1. For both deterministic or randomized algorithms, the error of simulating e ~ lHAt is no less than tt(At 3 ). 

Proof. Since deterministic algorithms are special cases of randomized algorithms, it is enough to consider randomized 
algorithms. Assume e ~ lHAt is simulated by U u with probability p w . Consider the Hamiltonians H\ and H 2 , in a 
given U w , let a±At, a 2 At, ■ ■ ■ , ax At be the total evolution time of Hi between two consecutive evolution of H 2 , while 
fliAt, (3 2 At, ■ ■ ■ ,/3xrAt are the total evolution time of H 2 between two consecutive evolution of Hi. For example, if 

jj^ _ e -iH 1 \ 1 At e -iH 2 \2At e -iH 3 \ 3 At e -iH 2 \4At e -iH 1 \ 5 At e -iH 4 \ 6 At e -iH 1 \ 7 At 

then ai = Ai, a 2 = A5 + A7 and (3\ = X 2 + A4. So, it is easy to see \K — K'\ < 1. Due to Lemma 1, the 
difference of trace distance is decided by E(\\ U u — Uq || 2 ) and || E(U U ) — Uo ||- If J2f=i a j 7^ 1 or J2f=i Pj 1 f° r 
some uj, || — Uo || 2 = £l(Ai 2 ), hence E(\\ — Uo || 2 ) = £l(Ai 2 ). Hence, we only need to consider the situation 
Y,f=i &j = 1 and Y^fLi Pj = 1 - Let us focus on the terms iHiH 2 H 1 At 3 and iH 2 H x H 2 At z . In e - iHAt , each of which 
has coefficients 1/6. If the simulation has an error less than 0(At 3 ), the coefficients of them in E{U {J f) must also be 
1/6, therefore the sum of them should be 1/3. However, we will show in every U u , the sum of these two coefficients 
is less than 1/3. Without loss of generality, assume K > K', let x 2 j-i = aj, for j — 1, • • • , K, and x 2 j — (ij, for 
j = 1, • • • ,K' . Then, the coefficient of iH\H 2 H\At z is the sum of XjXkXi, where j < k < I, j,l are odd, and k is 
even, while the coefficient of iHiH 2 HiAt 3 is the sum of XjXkXi, where j < k < I, j,l are even, and k is odd. Since 

^2fJ[ K Xj — J2f=i a j + SjLi Pj = 2, due to Lemma 2, the sum of the coefficients is always less than 1/3. Since 
there always exists 9(Ai 3 ) term in any simulation, the error of any simulation is no less than Q(At 3 ). □ 

Therefore, we have the main theorem. 
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Theorem 2. Any deterministic or randomized algorithm simulating e lHt by approximating it by Ylj=i e lAjtj , 
tj > 0, with error e satisfies 

N = n(t 3 / 2 e-^ 2 ). 

Proof. Assume the simulation is comprised of K stages, and in the j'-th stage, there are constant exponentials used 
to simulate e~ lHtj , where YljLi = From the above theorem, the final error is ^(X)j=i*|)j the minimum of 
which is Q(jfs). Hence, to grantee the final error is bounded by e, K = f^ 3/2 £~ 1/2 )- Therefore, N = tt(K) = 

n(< 3/2 £ -l/2)_ □ 

From Theorem [21 it is straightforward to obtain the following two Corollaries. 

Corollary 1. The deterministic algorithm based on the Baker-Campbell-Housdorf formula (Strang splitting) is asymp- 
totically optimal. 

Corollary 2. The randomized algorithm Algorithm 2 is asymptotically optimal. 



IV. SUMMARY 



In summary, we provide the randomized model of quantum simulation, and provide some randomized algorithms 
which are easier to implement than certain deterministic algorithms, but have the same efficiency. Moreover, we 
provide a lower bound for quantum simulation, therefore prove the optimality of the deterministic algorithm based on 
Strang splitting and one of our randomized algorithms. Note that, the lower bound and the optimality is under the 
assumption tj is positive in the simulation. Without this restriction, some algorithms have faster running time than 
the lower bound Furthermore, randomized algorithms also bring certain benefits in this unrestricted situation. 
For instance, when m = 2, to simulate e~ lHAt , it needs at least 7 exponentials to obtain an error bound 0(At 3 ) 0, 
however a randomized algorithm can obtain the same error bound with only 4 exponentials fl4j . 

We are grateful to Anargyros Papageorgiou, Joseph F. Traub, Columbia University for their very helpful discussions 
and comments. 
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